clear all
close all

[fname,fpath]=uigetfile('*.*','select the file to load'); %%% Get node data
VZDL=dlmread(strcat(fpath,fname),',',1,0);



badind=find(VZDL(:,7)==0);% removes "detections" for tag ID zero that mess up lambdas fit
VZDL(badind,:)=[];% removes "detections" for tag zero that mess up lambdas fit

%%%%%% Initialize Variables%%%%%%%%
% 2011 tag list %realtags=[10:1:59]; %This is the list of tags actually released, if this is not correct then this will not work!!!!
realtags=[10:1:29]; %2010 This is the list of tags actually released, if this is not correct then this will not work!!!!
month=VZDL(:,1);
day=VZDL(:,2);
year=VZDL(:,3);
hour=VZDL(:,4);
minute=VZDL(:,5);
seconds=VZDL(:,6);
tagID=VZDL(:,7);
tags=unique(tagID); % calls list of only tags ids heard
window=julianday(0,0,0,0,...
    str2double(inputdlg('Please enter window duration in minutes')),0); % Defines window length in minutes;
%samples=zeros(size(tags));
counts=zeros(size(tags));
correctedcounts=[];
dtxnTH=str2double(inputdlg('Please enter detection threashold. Ex: 0.9999'));
jday=zeros(size(day));
caphist=cell(size(tags));
lambdas=zeros(size(tags));
siteCode=str2double(inputdlg('Please enter the site code (1-7)')); % Site Codes: 1 = Pates, 2 = Cianbro, 3 = WWTP, 4 = PRSC, 
% 5 = GS,  6 = VZSRX400  7 = VZDL



for i= 1:length(day)
    jday(i)=julianday(year(i),month(i),day(i),hour(i),minute(i),seconds(i));
end
%set up the waitbar
wh=waitbar(0,'calculating noise levels')

%fore each entry in the detections array, count the number of times that
%tag occured within the window, and report the greatest count for each tag
%code in the "counts" array
for j=1:length(day)
    thistag=tagID(j);
    currenttime=jday(j);
    endtime=currenttime+window;
    ind=find(jday>=currenttime & jday <= endtime);
    thiscount=sum(tagID(ind,1)==thistag); %sum of occasions in window that tagID=thistag
    tagind=find(tags==thistag);
%     counts(tagind)=counts(tagind)+thiscount;
%     samples(tagind)=samples(tagind)+1;
 
        if thiscount>counts(tagind)
            counts(tagind)=thiscount;
        end
    waitbar(j/length(day),wh);
end
% %ind=find(samples==0);
% samples(ind)=1;
% counts=counts./samples;
close(wh);
bar(tags,counts);
bar(tags,counts);
title('Maximum number of occurences within window for each tag');
ylabel('Occurences');
xlabel('Tag Number');
%seperate the events count for tags that were released from tags that were
%not, so that a noise floor can be estimated using tags that are known to
%be 100% noise
noiseCounts=counts;
for k=1:length(realtags)
    realind=find(tags==realtags(k));
    if ~ isempty(realind)
        noiseCounts(realind)=0;%set the noise value to 0 for real tags
    end    
end
% 2011 take all noise, not just greater than 1- %noiseind=find(noiseCounts>0); %find the indicies of all noise tag codes
noiseind=find(noiseCounts>1); %2010 - find the indicies of all noise tag codes with noise greater than 1
noiseTags=tags(noiseind); %create a vector that just has the noise tag codes
noiseCounts=noiseCounts(noiseind); %create a vector that just has the counts from noise tags
averagenoise=mean(noiseCounts(3:length(noiseCounts)));%compute the average noise value for the whole signal, not including tags 1 and 2
lambdas=counts;        
lambdas(3:length(lambdas))=averagenoise; %use what figured out with counts to make array of lamdas for each receiver:will need to tweak later
%fitting the noise count to a negative exponatial equation
% lambda2 is a regression fit for negative exponential function


%this section of code removes counts of one, not sure why we did this,
%this is commented out now
% fitInds=find(noiseCounts~=1);
% fitTags=noiseTags(fitInds);
% fitCounts=noiseCounts(fitInds);
fitTags=noiseTags;
fitCounts=noiseCounts;
[b,bint,r,rint,stats]=regress(log(fitCounts),[ones(size(fitTags)),(fitTags)]);
lambdas2=exp(b(1)+b(2).*tags);
lambdas2(1:2)=counts(1:2);
if(max(lambdas2(3:length(lambdas2))>max(lambdas2(1:2))))
    meanNC=(mean(noiseCounts(3:length(noiseCounts))));
    lambdas2=ones(size(tags)).*meanNC;
    badInd=find(lambdas2<1);
    lambdas2(badInd)=1;
    lambdas2(1:2)=counts(1:2);
end
% begin section to make plots

figure;
subplot(1,2,1);hold on;
plot(tags,lambdas2,'r--');bar(tags,counts);
title('Counts, lambdas2 values, and cuttoff thresholds')

cutOffs=zeros(size(lambdas2)); % calculates cut off value for number of hits needed to NOT be random
for i=1:length(lambdas2)
    cutOffs(i)=poissinv(dtxnTH,lambdas2(i));
end
plot(tags,cutOffs,'m-')
legend('lambdas2','Counts','Threshold Value')
hold off


%now calculate lambda3, which is the bounded fit MAY WANT TO COMMENT OUT
%LATER
if(length(fitCounts)>10)
    logCounts=log(fitCounts);
    lowy=mean(logCounts(1:2));
    lowx=mean(fitTags(1:2));
    highy=mean(logCounts(length(logCounts)-5:length(logCounts)));
    highx=mean(fitTags(length(fitTags)-5:length(fitTags)));
    slope=(highy-lowy)/(highx-lowx);
    intcpt=lowy-(lowx*slope);
    lambdas3=exp(intcpt+slope.*tags); % lambdas 3 is an array of background noise values generated by a bounded exponential fit.
else
    lambdas3=ones(size(tags)).*str2double(inputdlg('enter minimum # detection for legit hit here for lambdas3')); % sets manual minimum # Detections in window for legtimate detection to 3 (user defined)
    lambdas3(1:2) = counts(1:2);
    
end
%begin plotting
subplot(1,2,2);hold on;
plot(tags,lambdas3,'r--');bar(tags,counts);
title('Counts, lambdas3 values, and cuttoff thresholds')
cutOffs=zeros(size(lambdas3));
for i=1:length(lambdas3)
    cutOffs(i)=poissinv(dtxnTH,lambdas3(i));
end
plot(tags,cutOffs,'m-')
legend('lambdas3','Counts','Threshold Value')
hold off



varString=inputdlg('Enter the variable you want to use to calculate the poisson threshold')% prompts which to use lambda 2 or 3

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wh=waitbar(0,'Marking Detections')
for j=1:length(day)
    thistag=tagID(j);
    currenttime=jday(j);
    endtime=currenttime+window;
    ind=find(jday>=currenttime & jday <= endtime);
    thiscount=sum(tagID(ind,1)==thistag); %sum of occasions in window that tagID=thistag
    tagind=find(tags==thistag);
    countThresh=poissinv(dtxnTH,eval(char(strcat(varString,'(tagind)')))); % returns minimum number of dtxns in window necessary to reach non random prob dtxnTH
   % 'eval' executes command contained in a string
    if thiscount>=countThresh 
        caphist{tagind}=[caphist{tagind},currenttime];
    end
    waitbar(j/length(day),wh);
end
close(wh);
[fname,fpath]=uiputfile('*.csv','select basepath for detections files');


for l=1:length(tags)
    if ~ isempty(caphist{l})
       capNum=length(caphist{l});
       tagVec=ones(capNum,1).*tags(l);
       siteVec=ones(capNum,1).*siteCode;
        dlmwrite(strcat(fpath,fname(1:length(fname)),'_detections_',num2str(tags(l)),'.csv'), [caphist{l}',tagVec,siteVec],'delimiter',',','precision',10);
    end
end

        
        
